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We study the behavior of two spatially distributed (sand- 
pile) models which are weakly linked with one another. Using 
a Monte-Carlo implementation of the renormalization group 
and algebraic methods, we describe how large-scale correla- 
tions emerge between the two systems, leading to synchro- 
nized behavior. 



I. INTRODUCTION 

Interacting systems have long been the subject of con- 
siderable interest and study, and they can display a rich 
variety of behaviors. One of the premier concepts that 
has emerged from such studies is the notion of synchro- 
nization, both in its traditional sense as well as in one 
of its modern variants (e.g., phase synchronization, lag 
synchronization, etc.) A particularly interesting 

class of systems to consider in light of these broaden- 
ing notions of synchronization is provided by automata 
("sandpile") models. Automata offer a rich assortment 
of well-studied, complex behaviors (e.g., self-organized 
criticality [|j), and have been used extensively in the lit- 
erature to model a variety of physical phenomena (e.g., 
1^). If two such automata are permitted to weakly in- 
teract 1^] , an interesting type of synchronization effect is 
seen to emerge: While small events on one sandpile are 
essentially uncorrelated with small events on the other, 
large-scale events are so highly correlated that not only 
is a large event on one sandpile almost always concomi- 
tant with a large event on the other sandpile, but the 
two events are in fact approximately equal in magnitude 
(with rms fractional deviation approaching zero). (See 
Fig. 1.) This result holds despite the weakness of the 
coupling between the sandpiles. Note that this "syn- 
chronization" between sandpiles is not periodic (i.e., the 
time interval between synchronized large-scale events is 
not fixed) , nor is it completely random (since correlations 
exist between temporal spacing of events and event size.) 

We can glean some basic insight into the origin of this 
form of synchronization from a relatively simple plausi- 
bility argument: As a large avalanche sweeps across one 
sandpile, it will, owing to the weak coupling, spill some 
small yet nontrivial number of grains onto the other sand- 
pile. Since sandpile models (like other self-organized crit- 
ical systems) are capable of generating avalanches of all 
sizes, these spilled grains could conceivably induce a large 
subsidiary avalanche in the second sandpile. In turn, this 
subsidiary avalanche could spill some grains back onto 
the first sandpile, and so on. In this manner, it is possible 
to imagine how high levels of correlations might develop 



between the two sandpiles during large events. It is thus 
reasonable to conjecture that through such feedback and 
mutual reinforcement, a large avalanche starting on one 
sandpile would have a very high probability of inducing 
a (simultaneous) avalanche of comparable magnitude on 
the other sandpile, despite the weakness of the coupling. 
Indeed, from this scenario one might also infer that per- 
haps these avalanches would not merely be comparable 
in size, but in fact nearly equal in size (that this should 
be the case is certainly plausible, though, admittedly, not 
compellingly obvious). 

While this intuitive argument is helpful, understanding 
the actual process by which inter-sandpile correlations 
develop proves to be surprisingly subtle and interest- 
ing. The purpose of this paper is to examine the nature 
of this complex interplay between interacting automata, 
and to show how it produces the observed large-scale syn- 
chronous behavior. To do so, we use a modification of a 
renormalization-group procedure originally developed by 
|p[-p^ for single-sandpile models, along with an algebraic 
technique. The renormalization procedure in fact turns 
out to be interesting in its own right, since it is imple- 
mented using a Monte Carlo method which proves to 
be highly efficient, thus rendering previously intractable 
renormalization calculations easily computable. As a re- 
sult, the methodology we employ here is likely to be ap- 
plicable to a variety of related automata problems. This 
paper is organized as follows. In Section || we present 
a prototype interacting-sandpile model and describe nu- 
merical simulations which demonstrate the emergence of 
large-scale synchrony. Section III contains a detailed dis- 
cussion of the renormalization procedure itself and its 
predictions. We also describe an alternative algebraic 
approach which proves useful for understanding certain 
key features of the model, including the appearance of 
so-called "coupling symmetry." Generalizations of the 
basic model are also described. 



II. BASIC MODEL AND PHENOMENOLOGY 

To begin, we recall a classic sandpile model studied 
by Dhar and Ramaswamy . The system consists of a 
two-dimensional square lattice, where to each lattice site 
ij one ascribes a non-negative integer h(i,j). The func- 
tion h{i,j) is called the "height" and represents the num- 
ber of "sandgrains" on a given site. The system evolves 
as follows: A lattice site is selected at random, and one 
grain is added to that site. Provided the new height does 
not exceed a certain critical value (taken througout this 
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paper to be 4), then nothing further happens. If, how- 
ever, the critical height is exceeded, then that site will 
"topple" {h{i,j) — > h{i,j) — 2) and spill one grain to its 
neighbor on the right (h{i,j + l) h{i, j + + and one 
to its neighbor above {h{i — l,j) h{i — l, j) + l). The af- 
fected sites on the right and above may in turn topple (if 
their heights are above the critical threshold), and so on. 
In this manner, it is possible for an avalanche to spread 
across the lattice. This type of model is referred to as 
a "directed" sandpile, since an avalanche can only prop- 
agate upwards or to the right. Once an avalanche has 
exhausted itself, a new grain is dropped onto a randomly 
selected site (from either sandpile), and the process re- 
peats. (We mention here that the asymptotic behavior 
of this model remains unchanged under a wide choice of 
dropping rules, as analyzed in Jill.) The dynamics of 
this model (which is representative of a large class of re- 
lated automaton models) is surprisingly complex and has 
been well-documented [|ll],0. One of its key features is 
that it exhibits avalanches of all sizes, where "size" refers 
to the total number of lattice sites that topple upon the 
addition of a single grain to the system. 

Now consider a system of two (independent) directed 
sandpile models, each evolving according to the rules out- 
lined above. The two sandpiles are assumed to be of iden- 
tical dimension, so that for every lattice site on the first, 
one can associate a corresponding lattice site on the sec- 
ond. (Visually, we imagine two planar lattices, one atop 
the other. Each site on the top sheet is matched with 
the site immediately below.) We'll refer to these two lat- 
tices as "Sheet A" and "Sheet B," respectively. We now 
allow the two sheets to interact according to the follow- 
ing rule: If a site on a given sheet topples, it will, as 
before, always spill two grains onto its own sheet (one 
to each of its neighboring sites on the right and above) . 
Now, however, we assume this toppling site also has some 
nonzero probability p of spilling an additional two grains 
onto the other sheet (one to each of the neighboring sites 
on the right and above). In this latter case we require 
h(i,j;s) — > h{i,j;s) — 4, where s denotes the sheet on 
which the toppling site lies, so that sand is "conserved." 
(Since this model is directed, it does not matter if the 
updating rules for the lattice are implemented sequen- 
tially or in parallel.) Note that when p = 0, the two 
sheets are dynamically independent. Our study focuses 
on the weak-coupling regime (p << 1). We'll henceforth 
refer to this particular model as the "two-sheet model" 
(generalizations will be described later). 

In a series of numerical simulations on this coupled 
system, we added grains (one at a time) to randomly 
selected sites (on either sheet), and monitored the re- 
sulting avalanches. (Most simulations were carried out 
for L = 1000 and very few avalanches in the simulations 
reached the edge, so we expect edge effects to be min- 
imal.) For each avalanche, we tracked the number of 
sites that toppled in each of the two sheets {Na,Nb), 
explicitly counting multiplicities if a given site under- 
went multiple topples. A representative graph is shown 



in Figure 2. The high density of datapoints along the 
X- and y- axes of Fig. 2 for small avalanche sizes indi- 
cates that these small avalanches remain largely confined 
to the sheet on which they started; only rarely will they 
spill over to the other sheet. This is not surprising, since 
the two sheets are only very weakly coupled to one an- 
other (p = 0.05) and thus the dynamics on each can 
be expected to be essentially independent. However, for 
large avalanches, a new trend is clearly seen to emerge: 
Even though, at each individual lattice site, the probabil- 
ity of a grain spilling over to the other sheet remains very 
low, nonetheless a large avalanche starting on one sheet 
is seen to have an equally large effect on the other sheet. 
In particular, the total number of sites that topple on 
each sheet (in a given avalanche) become nearly equal in 
magnitude (Fig. 2), with a root-mean-square fractional 
deviation that approaches zero (Fig. 3, solid line). Qual- 
itatively, it is as though the effective coupling strength 
between the two sheets increases with spatial scale. (We 
will return to this point later.) We note here (for empha- 
sis) that had we included in Fig. 2 only those avalanches 
that were touched off by the addition of grains to, say. 
Sheet A only, then the prominent trend towards the di- 
agonal seen in the figure for large avalanche size would 
be unaffected. 

Our goal is to show how this type of "large-scale syn- 
chrony" {LSS for short) arises between two such weakly 
interacting systems. We mention here that LSS also ap- 
pears in a larger class of models than just the numeri- 
cal example described above. For instance, one can con- 
struct "generalized two-sheet models," in which sites are 
allowed to spill either one or two grains onto neighboring 
sites on either/both sheets according to some probabil- 
ity matrix. As we'll show later, these generalized sys- 
tems (subject to some mild restrictions, namely, an over- 
all right/above symmetry) fall into the same universality 
class as our original two-sheet model and hence also ex- 
hibit LSS. (We have in fact found that power-law behav- 
ior - which is a characteristic of all the sandpile models 
to be discussed in this paper - is indeed not essential for 
the appearance of LSS); 

III. RENORMALIZATION-GROUP ANALYSIS 

The fundamental behavior of these weakly interacting 
automata can be understood using a renormalization- 
group analysis, as we now describe. We base our work on 
the renormalization procedure developed by Hasty and 
Wiesenfeld for (single sheet) directed-sandpile mod- 
els, which extended key work by Pietronero et al. [^,^. 
Adapting this procedure for systems of interacting au- 
tomata, we show how it can be used to explain the emer- 
gence of LSS. 

The basic idea behind the renormalization method is to 
repeatedly coarse grain the automaton lattice into cells 
of successively larger sizes. In our model, this means 
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grouping the lattice sites on each sheet into 2x2 blocks, 
then 4x4 blocks, then 8x8, and so on. At each stage, 
dynamical evolution rules are constructed which describe 
the behavior of the cells. Each time the basic cell size is 
increased, new dynamical evolution rules are constructed. 
The so-called "RG map" is a mapping that links the evo- 
lution rules for these different cell sizes. The behavior 
of the original automaton model on large spatial scales 
can then be deduced by examining the limiting behav- 
ior (i.e., fixed points) of this RG map. Before proceed- 
ing, we remark here upon an important distinction be- 
tween the course-graining procedure used here for our 
two-automata model and the procedure used for single- 
sheet models, as described by 0,D- Specifically, because 
we are interested in how each sheet behaves individu- 
ally, we do the coarse graining on each sheet individu- 
ally, rather than following the standard procedure which 
would naturally treat the full two-sheet lattice as a sin- 
gle entity and coarse grain it into cells which span both 
sheets. (If this latter procedure is followed, one finds that 
under the RG mapping, the two-sheet model converges 
to the single-sheet model of [ p^|jTl| | , proving that the two 
models have the same critical exponent 113]. Unfortu- 
nately, all information regarding correlations between the 
two sheets is lost.) 

In our model, the procedure works as follows. Imag- 
ine that the sites on Sheets A and B have already been 
coarse grained n times, so that the individual sheets are 
divided up into large "cells," each comprised of 2" x 2" 
individual lattice sites. We will use the term "cell pair" 
to collectively refer to a cell together with its associated 
cell on the other sheet. Adapting the renormalization 
scheme of ]^-p^] to our system, the evolution rules for a 
cell pair can be expressed in terms of a 3 x 3 probability 
matrix P". In particular, if a grain is added to a cell on a 
particular sheet, then the associated probability matrix 
is P'^ — ^ [a, 13 G {0, 1, 2}), where a is the number of 
grains spilled onto the same sheet (where the grain was 
added) and (3 is the number of grains spilled onto the 
other sheet. If a = 1 (or (3=1) then the direction of 
the spilled grain (up or right) is chosen at random. If 
a = 2 (or (3 = 2) then one grain is spilled up and the 
other to the right. For example, in our original two-sheet 
model, only spills of the types -Pg q, ^20 ^^^"^ ^22 '^^^ 
occur. (This characterization of a cell's dynamics is not 
quite complete. It is also necessary to distinguish be- 
tween two subcases of P^i, namely, the case when the 
spills on Sheets A and B are in the same direction (i.e., 
both to the right or both up), and when they are not 
(i.e., one to the right and one up). We'll denote these 
symmetric and anti-symmetric subcases as P^i g, Pii ^, 
respectively (Pi"i,, + P^^i,, = Pi"i).) 

The next step is to course grain the cells again (into 
2n-i-i ^ 2"+-'^ blocks), and construct the corresponding 
evolution rules governing the new, enlarged cell pairs. In 
other words, we wish to determine the RG map that re- 
lates P" to P""*"^. To do so, we utilize the procedure 



developed in ||lO|] for the case of a single-sheet automa- 
ton model. This method involves considering two-by-two 
blocks of the smaller (2" x 2") cells, and going through 
all possible combinatoric possibilities to derive the prob- 
abilities for the enlarged cells. We refer the reader to |Tc| ] 
for a description of this basic method. There is, however, 
one critical departure that we make from the procedure 
cited in . Namely, construction of the RG map for the 
single automaton case in Q proved arduous but analyti- 
cally tractable (the RG map contained on the order of 100 
terms). However, for our case of two weakly interacting 
automaton, the resulting RG map is much more complex 
(it contains several orders of magnitude more terms!), 
rendering its explicit calculation infeasible. As described 
below, to surmount this difficulty we use Monte Carlo 
simulations to numerically sample the various combina- 
toric possibilities associated with P" , and in this manner 
can approximate the probability matrix of the enlarged 
cells P"'^^. We then repeat this procedure and look at 
the limiting behavior of the resulting sequence of proba- 
bility matrices. 

Specifically, the renormalization mapping is computed 
as follows. Assume that the evolution matrix P", which 
describes the cell-pair dynamics for cells of size 2" x 2", 
is known. We now consider enlarged cells of size 2"+^ x 
2"+^, formed by grouping together four (2" x 2") cells 
into 2x2 blocks. The rules governing the behavior of 
these enlarged cells are obtained in the following man- 
ner. Imagine dropping a single grain onto the lower-left 
subcell of a (2 X 2) cell. For sake of argument assume 
this cell lies on the top sheet (Sheet A). The subcell will 
respond according to rules defined by the evolution ma- 
trix P". For example, the probability of that subcell not 
toppling is given by Pqq, while the probability of that 
subcell toppling onto all four of its downstream neigh- 
bors (two on each sheet) is given by P2 2- We continue 
to follow the avalanching process until all subcells (in 
both the cell on Sheet A and the corresponding cell on 
Sheet B) are quiescent. We now check where grains have 
exited the large cell. For example, if no grains left the 
large cell on Sheet A, even if there were some internal 
spills, then we have an event of type (0,0), while if two 
grains exited onto sheet B, one in each direction and none 
on sheet A, then we have an event of type (0, 2). How- 
ever, if two grains exited on sheet B, both in the same 
direction (say, to the right), and none on sheet A, then 
we have an event of type (0,1). Thus, we only count 
the number of directions in which grains exit the large 
cell, not the total number grains in each direction, for 
this part of the analysis. We then repeat this procedure 
a large number, of times. (Typically about 10^ trials 
are necessary for adequate accuracy.) At the end of this 
procedure we have a unnormalizcd matrix of evolution 
numbers, A^,/? which is the total number of type (a, /?) 
events which occurred. However, as discussed in [p|,p^, 
we need to "normalize" these probabilities in a specific 
manner. In order to do this we compute Aq^q differently 
than the other elements of the matrix A. The procedure 
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we use is that for each sample we take (i.e., for each drop 
onto the initial subcell), we count the number of grains 
that exit the large cell. If this number is larger than zero 
we add one less than this number to Ao,o- This "nor- 
malization" is very similar to the one used by Hasty and 
Wiesenfeld but slightly easier to compute in simulations 
(but would be more difficult analytically), and also easier 
to generalize for more complex sandpiles (such as ones in 
higher dimensions). In fact, when applied to the single 
sheet model studied by Hasty and Wiesenfeld, it leads to 
slightly more accurate estimates of the critical exponent 
than their procedure does. (The difference between the 
two procedures arises when an upper-right subcell spills 
two grains in the same direction out of the large cell. In 
this case our procedure adds one more to Ao,o than Hasty 
and Wiesenfeld's would.) 

Given the matrix A it is straightforward to compute 
Let |A| be the sum of all the elements of A. We 
view the elements as probability amplitudes and thus we 
need to convert them into true probabilities to continue 
the renormahzation procedure, thus, p = Aq_^/|A|. 

Representative results for the renormahzation process 
are as follows (accurate to about ±0.002): 
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This has three immediate consequences: 

(a) Having set the initial probability matrix P° to cor- 
respond with our original two-sheet model with weak 
coupling (p = 0.05), we find that the RG map quickly 
converges to the limiting probability matrix, P°° . The 
principal result here is that this limiting matrix is diago- 
nal. This shows that in any large avalanche (i.e., on large 
spatial scales) the number of topplings on the top sheet 
and bottom sheet will be approximately equal, thereby 
establishing the emergence of LSS in this model. 

(b) If we instead vary the values of the starting prob- 
ability matrix P" (corresponding to the generalized two- 
sheet models), we find that all (nontrivial) choices of 
starting configurations P'^ display the same limiting be- 
havior P°° in the renormahzation analysis. Hence these 
generalized models fall into the same universality class 
as the original, and therefore will all exhibit the same 
behavior (LSS) on large spatial scales. In particular, this 
shows, for example, that our original two-sheet model 
with weak coupling {p — 0.05 << 1) is in the same uni- 
versality class as a two-sheet model with full coupling 
{p = 1). In other words, when viewed on larger and 
larger spatial scales, the weakly interacting automata be- 
gin to act as though they were very strongly coupled. 
This strengthening of effective coupling constant with 



length scale can be regarded as the source of the high 
level of correlation between the two systems. 

c) If we examine the intermediate stages of the tran- 
sition process P'^ ^ ^ ...P°°, an interesting new 
feature appears: Under the RG map a general start- 
ing matrix P'^ will first become approximately symmetric 
(e.g., P^) prior to becoming nearly diagonal (e.g., P^^). 
(In the symmetric phase, P^^ w PJja^ ''^^'^ P^\ia ~ 0). 
Hence, the renormahzation analysis leads to a prediction 
that, in our original two-sheet automaton model, adding 
a grain to a site on, say. Sheet A, has an equal likelihood 
of inducing an (intermediate-size) avalanche on Sheet B 
as on Sheet A, even though the local dynamics dictate 
that small avalanches are much more likely to occur on 
the sheet to which the additional grain was added than 
on the other sheet. This is surprising in that we started 
with a model in which the coupling between neighbor- 
ing lattice sites was highly asymmetric (in the sense that 
each site is strongly coupled with its neighbors on its 
same sheet but only weakly coupled with its neighbors 
on the other sheet (i.e., p = 0.05)), and yet we are led 
to the conclusion that on larger length scales the effec- 
tive inter-sheet coupling becomes equal in strength to the 
intra-sheet coupling. A type of large scale 'coupling sym- 
metry' has thus emerged. This prediction was tested and 
borne out by numerical simulations of the automaton, as 
illustrated (by the dashed line) in Fig. 3. 

We can gain further insight into the nature of this sta- 
tistical synchrony and, in particular, the onset of this 
coupling symmetry, by forgoing the above renormahza- 
tion approach and instead utilizing an algebraic argu- 
ment based on work by Dhar for an analogous model 
(see also Zhang [Q). We will take our original two-sheet 
model and calculate the two-point correlation function 
C{x,y), which describes the expected number of top- 
plings at site y, due to the avalanche caused by adding 
a single grain to lattice site x. What we will prove is 
that if two sites x and y are sufficiently far apart, then a 
symmetry in the correlation function C'{x,y) ~ C{x,y) 
develops, where y denotes the site corresponding to y but 
on the opposite sheet. This calculation will show that 
adding a grain to a given site on one sheet will induce 
the same expected number of topplings on some distant 
site on its own sheet as it will on the corresponding (dis- 
tant) site on the other sheet - despite the weakness of 
the coupling between the two sheets. (In what follows, it 
will be convenient to let xl,xb denote the the neighbor- 
ing sites immediately to the left or below a site x, on the 
same sheet.) 

First, define a toppling matrix —A{x,y), which speci- 
fies the average number of grains that will spill directly 
from a site x to site y in the event that x topples. 
(Note that here we count only direct spillage between 
the two sites, not grains that might spill from a; to ?; by 
way of intermediate sites.) For our original two-sheet 
model, we have -A{y,y) = -2{1 + p); -A{yL,y) = 
l;-A{yB,y) = l;-A{yZ,y) = p;-A{y^,y) = p. All 
other components of A(a;, y) are zero. As in Dhar [|l^, it 
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is straightforward to show that the toppling matrix and 
correlation function obey the following general relation; 
J2z z)A{z, y) — 5x,y For our model, the only terms 
in the toppling matrix which contribute to the summa- 
tion are the four neighboring sites of y. Thus, the re- 
lationship reduces to C{x,yL) + C{x,yB) + pC{x,yL) + 
pC{x,yB) = 2(1 + p)C{x,y). We observe, however, that 
this relation is precisely the formula for the probability 
that a certain random walk starting at site x will hit site 
y. In this random walk, at every step the walker is equally 
likely to go up or right, and switches between sheets with 
probability p/{l + p). It follows then that the probabil- 
ity that a walker starting on one sheet will end up on 
that same sheet k steps later is (1 -|- [(1 — p)/(l -I- p)]'^)/2. 
Since this approaches 1/2 for large k, we conclude that 
C{x; y) « C{x; y) for x and y sufficiently far apart. (More 
precisely, the fractional difference between these correla- 
tions scales like [(1 — p) / {1 + p)]'' , where k is the distance 
between sites x and y in the "taxicab" metric, assuming 
of course that y is reachable from x, else both correlations 
are 0.) Hence, this demonstrates that on sufficiently large 
spatial scales, the intra-sheet and inter-sheet coupling be- 
come equal. 

IV. CONCLUSIONS 

In summary, we have examined (in the context of a 
few specific examples) the nature of the complex correla- 
tions arising between weakly interacting automata, and 
have used a Monte-Carlo implementation of a renormal- 
ization group analysis to understand the appearance of 
large-scale statistical synchrony in these systems. Since 
both our methods of analysis and the properties of SOC 
systems are extremely robust (e.g., the extension of the 
algebraic analysis to more general models is straightfor- 
ward), wc believe that the types of intcr-sandpile cor- 
relations found here will likely be a generic feature of 
other weakly coupled SOC systems. In fact, preliminary 
analysis siiggcsts that these properties even arise in some 
automata models which do not exhibit SOC, such as dis- 
sipative models. We note that our Monte Carlo approach 
for studying the RG map turns out to be remarkably effi- 
cient and may in fact provide the key to applying rcnor- 
malization to more complex automaton models. 

Lastly, we remark that the emergence of strong statis- 
tical correlations described here raises a number of in- 
teresting questions, including (i) Is there some univer- 
sal scaling law describing how the length scale at which 
strong correlations arise varies with the inter-sheet cou- 
pling strength; and (ii) Might it be possible to recast this 
phenomenon as a type of phase transition that occurs 
with increasing spatial scale? 

We would like to thank D. Dhar for helpful comments 
on the manuscript. 
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V. FIGURE CAPTIONS 



Figure 1. A representative time series. Shown is the 
total number of topples on each sandpile for each in a se- 
ries of avalanches. Note that large peaks occur simultane- 
ously and are approximately equal in magnitude, while 
smaller peaks are relatively uncorrelated in both time 
and size. (The dataset was generated from an automata 
model described in Section ||. Note that for illustrative 
purposes, we have added one to the avalanche sizes in or- 
der to avoid singularities associated with the logarithmic 
scaling in the plots.) 

Figure 2. Large-scale synchrony. The number of top- 
ples on each sheet during avalanches in the two-sheet 
model with coupling parameter p = 0.05 is shown. Ob- 
serve that for large avalanche sizes strong correlations de- 
velop, with A'^ and Nb becoming approximately equal. 

Figure 3. The root-mean-square fractional deviation 
frms (solid line) between Na and Nb vs. avalanche 
size {N = Na + Nb) for the data shown in Fig. 2. 
The decrease in frms with size indicates that, on large 
length scales, the two sheets behave as though they 
were strongly coupled. A related phenomenon, 'cou- 
pling symmetry' (see text), is illustrated by the dashed 
curve showing the average fractional deviation {fave — 
{{^A — Nb)/{Na + Nb))), where the average is com- 
puted over only those avalanches that were initiated by 
the addition of one grain to Sheet A. Observe that such 
avalanches, if small, remain primarily confined to Sheet 
A (as expected) , while large ones divide equally between 
the two sheets (since fave 0). 
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